<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of dempe1</title>
  <meta name="keywords" content="dempe1">
  <meta name="description" content="DEMPE1  Demonstrate parameter estimation on a simple 2nd order LTI system.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">examples</a> &gt; <a href="#">parameter_estimation</a> &gt; dempe1.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../../menu.html"><img alt="<" border="0" src="../../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\examples\parameter_estimation&nbsp;<img alt=">" border="0" src="../../../right.png"></a></td></tr></table>-->

<h1>dempe1
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>DEMPE1  Demonstrate parameter estimation on a simple 2nd order LTI system.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> DEMPE1  Demonstrate parameter estimation on a simple 2nd order LTI system.

   This is a simple demonstration of how to use the ReBEL toolkit for parameter estimation on
   a simple 2nd order LTI system.

   See also
   DEMSE1, DEMJE1, GSSM_LTI1

   Copyright (c) Oregon Health &amp; Science University (2006)

   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
   academic use only (see included license file) and can be obtained from
   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
   software should contact rebel@csee.ogi.edu for commercial licensing information.

   See LICENSE (which should be part of the main toolkit distribution) for more
   detail.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>	ADDRELPATH  Add a relative path which gets expanded into a absolute path</li><li><a href="../../.././ReBEL-0.2.7/core/cdkf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = cdkf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">cdkf</a>	CDKF  Central Difference Kalman Filter  (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>	CVECREP  Column vector replicate</li><li><a href="../../.././ReBEL-0.2.7/core/ekf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ekf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ekf</a>	EKF  Extended Kalman Filter</li><li><a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>	GENINFDS  Generate inference data structure from a generalized state space model and user defined inference parameters.</li><li><a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>	GENSYSNOISEDS  Generate process and observation noise data structures for a given InferenceDS data structure</li><li><a href="../../.././ReBEL-0.2.7/core/kf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = kf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">kf</a>	KF  Kalman Filter</li><li><a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>	REMRELPATH  Remove a relative path (which gets expanded into a absolute path)</li><li><a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>	SRCDKF  Square Root Central Difference Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/core/srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srukf</a>	SRUKF  Square Root Unscented Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/core/stringmatch.html" class="code" title="function match = stringmatch(string1,string2)">stringmatch</a>	STRINGMATCH  Returns match > 0 if string1 and string2 match (string2 can be a cell array of</li><li><a href="../../.././ReBEL-0.2.7/core/ukf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ukf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ukf</a>	UKF  Unscented Kalman Filter</li><li><a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_lti1.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_lti1</a>	GSSM_LTI1  Generalized state space model for simple LTI system</li></ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
</ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">% DEMPE1  Demonstrate parameter estimation on a simple 2nd order LTI system.</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%   This is a simple demonstration of how to use the ReBEL toolkit for parameter estimation on</span>
0004 <span class="comment">%   a simple 2nd order LTI system.</span>
0005 <span class="comment">%</span>
0006 <span class="comment">%   See also</span>
0007 <span class="comment">%   DEMSE1, DEMJE1, GSSM_LTI1</span>
0008 <span class="comment">%</span>
0009 <span class="comment">%   Copyright (c) Oregon Health &amp; Science University (2006)</span>
0010 <span class="comment">%</span>
0011 <span class="comment">%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for</span>
0012 <span class="comment">%   academic use only (see included license file) and can be obtained from</span>
0013 <span class="comment">%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the</span>
0014 <span class="comment">%   software should contact rebel@csee.ogi.edu for commercial licensing information.</span>
0015 <span class="comment">%</span>
0016 <span class="comment">%   See LICENSE (which should be part of the main toolkit distribution) for more</span>
0017 <span class="comment">%   detail.</span>
0018 
0019 <span class="comment">%=============================================================================================</span>
0020 
0021 clc;
0022 clear all;
0023 
0024 fprintf(<span class="string">'\nDEMPE1 : This demonstration shows how the ReBEL toolkit is used for parameter estimation\n'</span>);
0025 fprintf(<span class="string">'         on a 2nd order LTI system. The scalar observation is corrupted by additive white\n'</span>);
0026 fprintf(<span class="string">'         Gaussian noise.\n\n'</span>);
0027 
0028 
0029 <span class="comment">%--- General setup</span>
0030 
0031 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../gssm'</span>);         <span class="comment">% add relative search path to example GSSM files to MATLABPATH</span>
0032 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../data'</span>);         <span class="comment">% add relative search path to example data files to MATLABPATH</span>
0033 
0034 <span class="comment">%--- Initialise GSSM model from external system description script.</span>
0035 
0036 model = <a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_lti1.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_lti1</a>(<span class="string">'init'</span>);
0037 
0038 
0039 <span class="comment">%--- Generate some data</span>
0040 
0041 N  = 600;                                               <span class="comment">% number of datapoints</span>
0042 X  = zeros(model.statedim,N);                           <span class="comment">% state data buffer</span>
0043 y  = zeros(model.obsdim,N);                             <span class="comment">% observation data buffer</span>
0044 
0045 pnoise = model.pNoise.sample( model.pNoise, N);   <span class="comment">% generate process noise</span>
0046 onoise = model.oNoise.sample( model.oNoise, N);   <span class="comment">% generate observation noise</span>
0047 
0048 X(:,1) = [1 0]';                                          <span class="comment">% initial state</span>
0049 y(1)   = model.hfun( model, X(:,1), onoise(1), []); <span class="comment">% observation of initial state</span>
0050 <span class="keyword">for</span> j=2:N,
0051   X(:,j) = model.ffun( model, X(:,j-1), pnoise(j-1), []);
0052   y(j)   = model.hfun( model, X(:,j), onoise(j), []);
0053 <span class="keyword">end</span>
0054 
0055 
0056 <span class="comment">%--- Ask the user which inference algorithm to use</span>
0057 
0058 ftype = input(<span class="string">'Type of estimator [ kf / ekf / ukf / cdkf / srukf / srcdkf ] ? '</span>,<span class="string">'s'</span>);
0059 
0060 <span class="keyword">if</span> ~<a href="../../.././ReBEL-0.2.7/core/stringmatch.html" class="code" title="function match = stringmatch(string1,string2)">stringmatch</a>(ftype,{<span class="string">'kf'</span>,<span class="string">'ekf'</span>,<span class="string">'ukf'</span>,<span class="string">'cdkf'</span>,<span class="string">'srukf'</span>,<span class="string">'srcdkf'</span>})
0061     error(<span class="string">'Unknown estimator!'</span>);
0062 <span class="keyword">end</span>
0063 
0064 
0065 <span class="comment">%--- Setup argument data structure which serves as input to</span>
0066 <span class="comment">%--- the 'geninfds' function. This function generates the InferenceDS and</span>
0067 <span class="comment">%--- SystemNoiseDS data structures which are needed by all inference algorithms</span>
0068 <span class="comment">%--- in the PiLab toolkit.</span>
0069 
0070 Arg.type = <span class="string">'parameter'</span>;                              <span class="comment">% inference type (parameter estimation)</span>
0071 Arg.tag = <span class="string">'Parameter estimation for GSSM_LTI1 system.'</span>;  <span class="comment">% arbitrary ID tag</span>
0072 Arg.model = model;                                   <span class="comment">% GSSM data structure of external system</span>
0073 Arg.paramFunSelect = <span class="string">'ffun'</span>;                         <span class="comment">% Only use gssm_lti's FFUN for parameter estimation purposes</span>
0074 Arg.paramFFunOutIdxVec = [1];                        <span class="comment">% Only use the first element of the output vector of FFUN as the</span>
0075                                                      <span class="comment">%  observation. This is for auto-regressive state space models.</span>
0076 
0077 InfDS = <a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>(Arg);                                <span class="comment">% Create inference data structure</span>
0078 [pNoise, oNoise, InfDS] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS,ftype); <span class="comment">% Generate process and observation noise data structures</span>
0079 
0080 <span class="comment">%--- Setup runtime buffers</span>
0081 
0082 Wh = zeros(InfDS.statedim,N-1);                      <span class="comment">% meta system state estimation buffer</span>
0083 Wh(:,1) = 0.01*randn(InfDS.statedim,1);              <span class="comment">% initial condition of state E[Wh(0)]</span>
0084 Pw = eye(InfDS.statedim);                            <span class="comment">% initial state covariance</span>
0085 U2 = X(:,1:end-1);                                   <span class="comment">% exogenous input for observation function of meta system</span>
0086 Y  = X(InfDS.paramFFunOutIdxVec,2:end);              <span class="comment">% meta system observations</span>
0087 
0088 
0089 <span class="comment">%--- Call inference algorithm / estimator</span>
0090 
0091 <span class="keyword">switch</span> ftype
0092 
0093     <span class="comment">%------------------- Linear Kalman Filter --------------------------------------</span>
0094     <span class="keyword">case</span> <span class="string">'kf'</span>
0095 
0096         pNoise.adaptMethod = <span class="string">'anneal'</span>;             <span class="comment">% process noise adaptation method : annealing</span>
0097         pNoise.adaptParams = [0.9 1e-8];           <span class="comment">% annealing parameters [anneal_factor minimum_variance]</span>
0098 
0099         [Wh, Pw, pNoise] = <a href="../../.././ReBEL-0.2.7/core/kf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = kf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">kf</a>(Wh(:,1), Pw, pNoise, oNoise, Y, [], U2, InfDS);
0100 
0101     <span class="comment">%------------------- Extended Kalman Filter ------------------------------------</span>
0102     <span class="keyword">case</span> <span class="string">'ekf'</span>
0103 
0104         pNoise.adaptMethod = <span class="string">'anneal'</span>;
0105         pNoise.adaptParams = [0.9 1e-8];
0106 
0107         [Wh, Pw, pNoise] = <a href="../../.././ReBEL-0.2.7/core/ekf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ekf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ekf</a>(Wh(:,1), Pw, pNoise, oNoise, Y, [], U2, InfDS);
0108 
0109 
0110     <span class="comment">%------------------- Square Root Unscented Kalman Filter -----------------------------------</span>
0111     <span class="keyword">case</span> <span class="string">'ukf'</span>
0112 
0113         alpha = 1;        <span class="comment">% scale factor (UKF parameter)</span>
0114         beta  = 3;        <span class="comment">% optimal setting for Gaussian priors (UKF parameter)</span>
0115         kappa = 0;        <span class="comment">% optimal for state dimension=2 (UKF parameter)</span>
0116 
0117         InfDS.spkfParams = [alpha beta kappa]; <span class="comment">% UKF parameter vector</span>
0118 
0119         pNoise.adaptMethod = <span class="string">'anneal'</span>;
0120         pNoise.adaptParams = [0.9 1e-8];
0121 
0122         [Wh, Pw, pNoise] = <a href="../../.././ReBEL-0.2.7/core/ukf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ukf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ukf</a>(Wh(:,1), Pw, pNoise, oNoise, Y, [], U2, InfDS);
0123 
0124 
0125     <span class="comment">%------------------- Square Root Central Difference Kalman Filter ---------------------------</span>
0126     <span class="keyword">case</span> <span class="string">'cdkf'</span>
0127 
0128         InfDS.spkfParams = sqrt(3);     <span class="comment">% scale factor (CDKF parameter h)</span>
0129 
0130         pNoise.adaptMethod = <span class="string">'anneal'</span>;
0131         pNoise.adaptParams = [0.9 1e-8];
0132 
0133 
0134         [Wh, Pw, pNoise] = <a href="../../.././ReBEL-0.2.7/core/cdkf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = cdkf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">cdkf</a>(Wh(:,1), Pw, pNoise, oNoise, Y, [], U2, InfDS);
0135 
0136 
0137     <span class="comment">%------------------- Square Root Unscented Kalman Filter -----------------------------------</span>
0138     <span class="keyword">case</span> <span class="string">'srukf'</span>
0139 
0140         alpha = 1;        <span class="comment">% scale factor (UKF parameter)</span>
0141         beta  = 3;        <span class="comment">% optimal setting for Gaussian priors (UKF parameter)</span>
0142         kappa = 0;        <span class="comment">% optimal for state dimension=2 (UKF parameter)</span>
0143 
0144         InfDS.spkfParams = [alpha beta kappa]; <span class="comment">% UKF parameter vector</span>
0145 
0146         pNoise.adaptMethod = <span class="string">'anneal'</span>;
0147         pNoise.adaptParams = [0.9 1e-8];
0148 
0149         Sw = chol(Pw)';
0150 
0151         [Wh, Sw, pNoise] = <a href="../../.././ReBEL-0.2.7/core/srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srukf</a>(Wh(:,1), Sw, pNoise, oNoise, Y, [], U2, InfDS);
0152 
0153 
0154     <span class="comment">%------------------- Square Root Central Difference Kalman Filter ---------------------------</span>
0155     <span class="keyword">case</span> <span class="string">'srcdkf'</span>
0156 
0157         InfDS.spkfParams = sqrt(3);     <span class="comment">% scale factor (CDKF parameter h)</span>
0158 
0159         pNoise.adaptMethod = <span class="string">'anneal'</span>;
0160         pNoise.adaptParams = [0.9 1e-8];
0161 
0162         Sw = chol(Pw)';
0163 
0164         [Wh, Sw, pNoise] = <a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>(Wh(:,1), Sw, pNoise, oNoise, Y, [], U2, InfDS);
0165 
0166     <span class="comment">%---</span>
0167     <span class="keyword">otherwise</span> error(<span class="string">'That estimator/filter type is not recognized.'</span>);
0168 
0169 <span class="keyword">end</span>
0170 
0171 
0172 <span class="comment">%--- Display results</span>
0173 
0174 true_model_trace = <a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(model.params,N-1);
0175 
0176 figure(1); clf;
0177 p11=plot(true_model_trace(1,:),<span class="string">'b--'</span>,<span class="string">'linewidth'</span>,3); hold on;
0178 p12=plot(true_model_trace(2,:),<span class="string">'r--'</span>,<span class="string">'linewidth'</span>,3);
0179 p21=plot(Wh(1,:),<span class="string">'b'</span>);
0180 p22=plot(Wh(2,:),<span class="string">'r'</span>); hold off;
0181 axis([1 N -1.5 2.5]);
0182 legend([p11 p21 p12 p22],<span class="string">'true parameter 1'</span>,<span class="string">'estimated parameter 1'</span>,<span class="string">'true parameter 2'</span>,<span class="string">'estimated parameter 2'</span>,0);
0183 ylabel(<span class="string">'model parameters'</span>)
0184 xlabel(<span class="string">'time'</span>);
0185 title([<span class="string">'DEMPE1 : LTI System Parameter Estimation  ( '</span> ftype <span class="string">' )'</span>]);
0186 
0187 
0188 <span class="comment">%--- House keeping</span>
0189 
0190 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../gssm'</span>);       <span class="comment">% remove relative search path to example GSSM files from MATLABPATH</span>
0191 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../data'</span>);       <span class="comment">% remove relative search path to example data files from MATLABPATH</span></pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>